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Abstract 

We study the homogeneous nucleation process in Stillinger- Weber silicon in the NVT ensemble. 
A clear first-order transition from the liquid to crystal phase is observed thermodynamically with 
kinetic and structural evidence of the transformation. At 0.75 T m , the critical cluster size is 
about 175 atoms. The lifetime distribution of clusters as a function of the maximum size their 
reach follows an inverse gaussian distribution as was predicted recently from the classical theory of 
nucleation (CNT). However, while there is a qualitative agreement with the CNT, the free energy 
curve obtained from the simulations differs significantly from the theoretical predictions, suggesting 
that the low-density liquid phase found recently could play a role in the nucleation process. 

PACS numbers: 64.70.Dv 64.60. Qb 82.60.Nh 
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I. INTRODUCTION 



The classical nucleation theory (CNT) has been extensively tested in systems with rel- 
atively simple two-body interactions such as colloids or globular proteins 1 i 2 i 3 i 4 . These 
molecules are large and move slowly, making it possible to follow the crystallization pro- 
cess experimentally using various techniques of microscopy. Moreover, these systems can 
also be represented accurately by theoretical models of hard- and soft-spheres, which can 
crystallize on numerical time scales. It is therefore possible to characterize fully the micro- 
scopic mechanisms responsible for nucleation in terms of the CNT, which works particular 
well for these systems. 

There has also been a number of studies going beyond the soft-sphere models. In par- 
ticular, there has been considerable work devoted to the nucleation of Lennard- Jones mod- 
e b 5 i 6 i 7 i 8 i 9 . Very little work has been done, however, on more complex materials such as 
oriented liquids — water or tetrahedral semiconductors, for example. Recently, Matsumoto 
et al. — , using considerable computing power, managed to follow one occurrence of crys- 
tallization in a 300 ns run of a 512-molecule simulation of water in the canonical ensemble 
at 230 K. Clearly more simulations are needed in water but also in simpler oriented liquids 
such as silicon, which shows a similar phases diagram around melting as both liquids show 
a temperature of density maximum and their density falls off by ~ 10% from the disordered 
liquid to the tetrahedral crystalline structure. As with water, there has been very few works 
studying nucleation in this technologically important material. 

Depending on the cooling rate, previous numerical work has shown that supercooled 
liquid silicon transforms in a glassy 12 ' 13 or amorphous 14 > 15 state. Recently, it was indicated 
that this transition takes place just below a liquid-liquid transitio n 16 ! 17 : at zero pressure in 
the Stillinger- Weber silicon, the low density liquid (LDL), which is thermodynamically and 
structurally contiguous to the amorphous solid, crystallizes rapidly (around 10 ns) at 1050 
K 16 i 17 whereas the more common high-density liquid (HDL) does not at any temperature on 
a simulation timescale. In order to circumvent the difficulty to crystallize /-Si, Uttormark 
and colleagues — embedded a spherical crystal seed containing 400-800 atoms in bulk liquid 
and analyzed the growth and dissolution of clusters. They found that the critical size for a 
crystallite to grow to macroscopic size was of 140 and 1400 atoms at 60% and 85% of the 
melting temperature (T m ). Working with a similar method, Bording and Tafto^ inserted 
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a crystallite in an amorphous matrix of 4096 germanium atoms and estimated the critical 
cluster radius to be 2 nm (around 1500 atoms) at 60% T m . 

In this paper, we show that liquid silicon can crystallize in the NVT ensemble on timescale 
accessible by MD simulation without going through the low-density liquid phase. We also 
show that the nucleation process, while qualitatively consistent with CNT, differs quantita- 
tively from it. 

The organization of the paper is as follow. We show the behavior of the thermodynamic, 
kinetic and structural properties during the phase transition in section IIII Al In section 
IIII B 1| we analyze the nucleation and crystallization process through the evolution of the 
cluster that will eventually crystallize the whole system in relation to CNT. Then, in section 
1111 B 21 we compute and compare the free energy of clusters between CNT and the simulation 
data. Finally, we look at the lifetime of small clusters in the supercooled liquid before 
nucleation takes place in section IIII B 31 

II. METHODS 

The molecular dynamical simulations (MD) for this work are performed in the canonical 
(NVT) ensemble at the K crystalline density, i.e. 2.32 g/cm 3 , and in the isothermal-isobaric 
(NPT) ensemble at pressure. All simulations are done at 1250 K (75% T m ) in a cubic box 
containing 10648 atoms, with periodic-boundary conditions. This size is sufficiently large 
to avoid catastrophic crystal growth due to interactions between the images of the critical 
crystallite, which is estimated to be around 200 atoms (see below). 

We use the extended-system method of Andersen to control pressur o 20 i 21 i 22 and Hoover's 
constraint method for the temperatur o 23 i 24 i 25 . Newton's equations of motion are integrated 
with a fifth-order Gear predictor-corrector and a time step At = 1.15 fs. Simulations are 
typically equilibrated over 50 000 At (58 ps) and data are accumulated over 10 6 At (several 
ns). Atomic interactions are represented by the Stillinger- Weber potential (SW), developed 
to reproduce accurately the crystalline and liquid state of Si 2 ^. 

Starting with a liquid well equilibrated at 2900 K, we generate nine independent trajec- 
tories in NVT conditions at 2.32 g/cm 3 and 75% of T m , a degree of undercooling similar to 
that used for a wide variety of materials both experimentally 2 ^ and numerically— Of these 
nine trajectories, six crystallize within 10 ns and are numbered 1 to 6; the fastest, simulation 
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# 1, crystallizes within 1.5 ns. 

Following previous work on liquid S i 11 ! 17 , we use as order parameter the smallest three- 
dimensional closed-ring structures that can be associated with a given crystalline lattice. 
These clusters, shown in Fig. are the smallest elementary building blocks for wurtzite, 
diamond and /3-tin structures and are defined topologically: the wurtzite lattice is associated 
with a 12-atom cluster composed of two sixfold rings connected at three points while the 
diamond and diamond can be described topologically by a single 10-atom cluster with four 
sixfold rings back to back. To establish the connectivity of these clusters, the first-neighbor 
cut-off is set to 2.75 A, a value similar to that used in these high-quality amorphous networks. 
This is somewhat shorter than the typically nearest-neighbor distance used in liquids (which 
is about 3.0 A) as it focuses on local crystalline order. 

These elementary clusters are present with a low density in the liquid (p C r y stai ~ 5 — 10 
% at.) as well as in high-quality amorphous models prepared using the modified WWW 
bond-switching algorithm (p crys tai ~ 1 — 5 % at.)2£. These blocks provide therefore a much 
more sensitive measure of crystallinity than the structure factor or the RDF. 

Our criteria are different from those used in a previous study of the nucleation of crystal- 
lites implanted into a SW liquid by Uttormark, Thompson and Clancy^. In this case, the 
description of a crystallite nucleus is defined uniquely based on a mixture of energetic, topo- 
logical and geometric constraints. For an atom to be part of a crystallite: (i) its three-body 
energy in SW potential of fourfold or fivefold coordinated atoms (within a 3.35 A nearest- 
neighbor distance) must be lower than 0.4336 eV; (ii) it must posses four nearest neighbors 
and at least three of them are also fourfold coordinated; (iii) its angular bond angles meet 
this criterion: 5^ =1 (cos @j + 1/3) 2 < 0.4 (where 0j is the angle between nearest neighbors of 
a fourfold coordinated atom.) The crystallites identified with this method are less compact 
than those flagged with our topological order parameter. This is particularly true for small 
crystallites (less than 20 atoms), which tend to be open and stringy, like twisted polymers, 
with Uttormark's criteria. The two methods converge, however, for larger clusters, near and 
beyond the critical size, where a clear definition of surface is less important. 
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FIG. 1: (color online) The three basic building blocks associated with the crystalline order param- 
eter. The wurtzite basic block (left) is a 12-atom cluster composed of two sixfold rings whereas the 
diamond basic block (middle) is a 10-atom cluster with four sixfold rings. The /3-tin basic block 
(right) is equal to a diamond basic block where the tetrahedra are compressed in one direction and 
elongated along the two others axes. 

III. RESULTS AND DISCUSSION 
A. Phase transition 

Homogeneous nucleation is often difficult to obtain numerically, especially in oriented 
solids such as Si and water which display a crystalline structure far from that of the liquid 
phase. It took months of computer time to simulate homogeneous nucleation in TIP3P 
water. Studies using SW Si failed to find traces of nucleation in a 5000-atom cell after a 
1-ns simulation.— 

In view of these results, and because classical nucleation theory (CNT) predicts that 
nucleation and crystallization is obtained more rapidly for strong undercooling and larger 
system size, we choose to simulate a larger cell, with more than 10 4 atoms, simulated over 
10 ns at 0.75 T m . 

As shown in Figure El this is sufficient to observe homogeneous nucleation, from the pure 
liquid phase, in the NVT ensemble. While the data presented in this figure are for simulation 
#1, a run that crystallizes particularly quickly, the overall properties of the transition are 
identical to run #2 to 6. The top curve shows a brutal drop in the potential energy of the 
system, from —3.75 to —3.95 eV/atom, indicating a clear thermodynamical transition after 
1.5 ns of simulation. The phase transition is also visible by following the change in pressure 
(bottom panel). As the density is maintained at the crystalline value, the pressure in the 
liquid phase is negative; it changes sign at the liquid- crystal transition since the crystal 
density at 1250 K and GPa is slightly lower than at K, and since the new structure 
contains grain boundaries. The liquid to solid phase transition is clearly seen in the kinetics 
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FIG. 2: Evolution of the energy (top), mean square displacement (middle) and pressure (bottom) 
during the liquid-crystal phase transition of Si with NVT conditions at 1250 K and 2.32 g/cm 3 . 
These results are for simulation # 1, which crystallizes the fastest. While the other simulations 
take longer to crystallize, their evolution is similar. 

of the system (middle panel): in the supercooled liquid, the diffusion is significant, with 
D = 5.4 • 10~ 6 cm 2 /s; it drops suddenly at the transition to become vanishingly small, a 
clear indication of a liquid to solid transition. 

Under the NVT conditions described in the introduction, the mean pressure of the su- 
percooled liquid is -1.9 GPa. In a previous worki£, we studied the transition from high 
density liquid (HDL) to low density liquid (LDL) in Stillinger- Weber Si and showed that 
this transition does occur at around 1250 K and -2 GPa but moves to lower temperatures 
as the pressure is increased. The current simulations are therefore slightly above the HDL 
to LDL transition, and we seem to observe a pure liquid-crystalline transition: the liquid 
before the transition has a RDF and a diffusion constant characteristic of the HDL and there 
is no trace of a LDL phase during the crystallization process. 
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FIG. 3: Radial distribution function before and after the liquid-crystal phase transition of Si in 
NVT conditions at 1250 K and 2.32 g/cm 3 . The RDF is characteristic of a crystalline state after 
the transition (2.7 ns) and of a liquid before the transition (0 ns). These results are taken from 
Simulation #1. 

Changes in the structural properties of this model as the transition occurs are shown 
in the next two figures. At t = ns, the radial distribution function (RDF) (see Fig. EJ) 
is typical of that of a liquid, with little structure beyond the broad second-neighbor peak. 
The nature of the RDF is totally different after the transition, with well-defined crystalline 
peaks up to 9 A and beyond. In the liquid phase, the system contains very few crystalline 
building blocks and p C r y stai fluctuates between 5 and 10 % of all the atoms (see Fig. HJ top 
panel). After the transition, more than 85 % of atoms belongs to a diamond and/or wurtzite 
crystalline blocks, with a probability higher for diamond structures except in trajectory 

The co-existence of two crystalline structures is not surprising since, with a cutoff of 
3.77 A, the SW potential cannot differenciate between the diamond and wurtzite structures 
at zero temperature: these two structures start to differ only at their third-neighbor shell, 
at 4.50 and 3.91 A, respectively. It is therefore only the thermal vibrations, bringing the 
third-neighbor shell atoms inside the cut-off from time to time, that allow the potential to 
distinguish between these two crystalline structures. With long enough annealing, we expect 
the wurtzite structures to disappear completely. 

For its part, the liquid phase is characterized by a low density of crystalline building 
blocks. Moreover, these crystallites tend to be small, counting less than 20 atoms, on 
average. Before crystallization begins, the number of independent nuclei oscillates between 
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FIG. 4: (color online) Proportion of atoms in elementary blocks (top) and number of independent 
clusters (bottom) during the liquid-crystal phase transition of Si in NVT conditions at 1250 K and 
2.32 g/cm 3 . The proportion of atoms in diamond and/or wurtzite crystalline structures (0 crystal) 
increases rapidly reaching a value close to 1 after the transition. These results are for Simulation 
# 1- 

40 and 50. As crystallization occurs, however, the largest nucleus grows rapidly, absorbing 
the smaller crystallites and forming a single system-size cluster; the number of independent 
crystallites decreases constantly during this process (FiglH bottom panel). 

B. Characterization 

1. Stability of crystallites 

It is possible to characterize more finely the crystallization by following the crystalline 
precursor as it takes over the simulation cell. This is achieved by following the evolution of 
all crystallites by steps of 1.15 ps. During this short time some crystallites appear, other 
vanish, while the rest might evolve significantly; a set of rules must therefore be established 
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TABLE I: Characteristic times of the crystalline precursor that gives rise to crystallization of the 
supercooled liquid. 



Simulations 


Time (ns) 


in 


^ZUU 


^IlUC 




^crys 


1 


0.06 


0.38 


0.46 


0.59 


1.60 


2 


0.64 


0.75 


0.84 


1.02 


2.75 


3 


3.14 


3.30 


3.32 


3.62 


4.85 


4 


3.30 


3.67 


4.20 


4.68 


6.00 


5 


5.14 


5.35 


5.32 


5.76 


8.00 


6 


7.79 


7.98 


7.98 


8.27 


9.75 



to identify uniquely and reversibly each aggregate: (1) At least three atoms must remain 
together over one time interval for a cluster to survive; a failing test indicates that the 
aggregate has dissolved. (2) When two or more crystallites merge together, the one with the 
highest number of surviving atoms is considered the progeny, the other one ceases to exist. 
(3) If, on the other hand, a cluster splits into multiple parts, the new aggregate containing the 
highest number of original atoms becomes the progeny and the other clusters are considered 
newborn. Using this analysis, we can then follow the evolution of the crystalline precursor 
by tracing back its ancestors. 

In order to compare between the 6 runs that crystallize, we separate the time evolution 
into four periods. The instant of birth of the crystalline precursor is defined as to (see Tab.QJ). 
From this time, it may take several hundreds of picoseconds (about 200 to 900 ps) for this 
embryo to reach a critical size, at time t nnc . The nucleation time t nuc is defined as the point 
in time where the size of the aggregate starts growing steadily, as seen in Fig. At this 
point, the system leaves the incubation regime to enter the steady-stade of nucleation and 
crystallization as such takes place. 

CNT predicts that a cluster of over-critical size should grow continuously whereas under- 
critical size crystallites tend to dissolve, in both cases, to lower their free energy. Statistical 
fluctuations can foil those predictions around the critical size, however, and move from 
under-critical to over-critical size and vice- versa. This explains why we define t nuc not as the 
first time when the cluster reaches the critical size, but the first time it reaches it for good. 
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FIG. 5: Evolution of the crystalline precursor during the liquid-crystal transition of Si in NVT 
conditions at 1250 K and 2.32 g/cm 3 . These results are taken from the first simulation (#1) to 
crystallize. 

For example, while at i nuc crystalline precursors are composed in average of 160 atoms, they 
have often reached a size of 200 atoms or more before. This fine characterization of t nuc is 
probably not needed, however. Looking at Table HI t nuc appears closely correlated with t 2 oo, 
the point in time where the crystallite reaches a maximum size of 200 atoms for the first 
time. The number of clusters reaching a 200-atom size or more and then dissolving into the 
liquid is extremely small. Thus, the critical cluster size should be around 175 atoms for Si 
at 1250 K, in agreement with the estimate of Uttormark et a/.—. 

^From t nuc , the crystallization per se proceeds rapidly into a steady growth regime which 
lasts about 2 ns. The crystallization time, t crys , is defined as the moment when the size of 
the largest cluster stops growing. 

For all simulations, it is possible to trace back the critical cluster to its appearance as a 
small aggregate of about 20 atoms, at to- By selection, this cluster should live longer than 
most other under-critical crystallites. As shown in Fig. the size of this cluster typically 
oscillates for a long time, aggregating and loosing atoms until it reaches a critical size at 
t nuc and then starts growing for good. 

Surprisingly, while the cluster size oscillates, its composition changes considerably. 
Throughout the incubation regime, the crystalline precursor changes its composition sig- 
nificantly: very few atoms of the original cluster remain part of it until the nucleation phase 
starts. In half the simulations, less than 50% of the original atoms are part of the cluster for 
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TABLE II: Proportion of atoms participating into the crystalline precursor permanently and 90 % 
of the time during the incubation and steady-state regime of nucleation. Starting from atoms who 
belong originally to the crystallite at time to until t nnc in the incubation phase and from t nnc until 
£500 in steady-state. The interval between each configurations snapshot is 1.15 ps. 



Simulations 


Persistence of atoms part of the crystalline precursor 




From t to t nuc 


From t nuc to t 500 




permanent 


90% of the time 


permanent 


90% of the time 


1 


18% 


82% 


52% 


81% 


2 


4% 


64% 


43% 


77% 


3 


38% 


79% 


18% 


57% 


4 


0% 


40% 


18% 


51% 


5 


0% 


46% 


13% 


40% 


6 


0% 


20% 


29% 


49% 



90% of the time interval between i and t nuc (see tab.|Hl- Even in the steady-growth regime, 
starting at i nuc , the crystallite continues to exchange atoms with the liquid. For most of the 
runs, less than half the 160 or so atoms present at t nnc remain in the clusters for 90% of the 
time in this interval until £500; as the growth takes place a significant fraction of the atoms 
move back and forth between the crystallite and surrounding. These results are in line with 
a previous study on growth and dissolution of implanted LJ crystallites with a critical size 
similar to that of our system 6 which shows that the probability of dissolution, while decrease 
rapidly with cluster size above the critical size is still non-negligible for clusters 50% bigger 
than critical size. Although we follow a cluster that will not dissolve totally, the considerable 
atomic exchange is a reflection of this tendency. While the critical aggregate's composition 
changes rapidly, its position remains almost fixed in space, its center of mass hardly moving 
except by aggregation. The crystalline precursor is therefore not a static crystalline seed 
slowly growing throughout the nucleation process; there is a constant exchange of matter 
with the surrounding liquid even for post-critical sizes. 
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FIG. 6: Evolution of nucleation and crystallization during the liquid-crystal phase transition of 
SW Si at 1250 K and 2.32 g/cm 3 . The configurations show atoms who belong only to crystalline 
structures at 0, 0.58, 0.86, 1.15, 1.44 and 1.73 ns respectively for simulation #1. 
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FIG. 7: Free energy (divided by fc#T) of crystallites as a function of their size in the NVT ensemble. 
The simulation data are computed from the equilibrium probability of presence for clusters with 
the basic blocks analysis (□) or the criteria of Uttormark et al. ^ (+). The CNT curve computed 
with the AF s i value from thermodynamic integration is indicated as Fit A. A better fit is given by 
Fit B. Details are discussed in the text. 

2. Free energy 

It is formally straightforward to compare the simulations with the predictions of CNT on 
the thermodynamics of crystal growth. The free energy curve of crystallites can be obtained 
from the simulations by plotting the equilibrium probability P eq (n) to find a crystallite of 
size n in the metastable liquid^^. 

We compute P eq (n) in the supercooled liquid, accumulating data until the largest cluster 
reaches 500 atoms, less than 5 % of the total number of atoms but over the critical size (see 
Tab. H}, and over all runs. This distribution is directly connected to the free energy AF(n) 
associated with these clusters: 



where N(n) is the number of clusters of size n present in the liquid, ksT is the Boltzmann 
constant times temperature and C, a constant. 

The CNT offers another way to compute the free energy. In a simple relation, the energy 
gain in the formation of a new phase is balanced by the cost to produce an interface between 
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the old and new phases: 

AF(n) = AF sl -n + a- n 2/3 (3) 

where AF S ; = F s — Fi is the Helmholtz free energy difference between solid and liquid states 
in NVT conditions, a = A ■ 7 with 7 the surface tension, A = (3Qtt / pl) 1 ^ for spherical 
crystallites and p s the density of the solid phase. While the Helmholtz free energy difference 
AF s i is relatively easy to obtain, the evaluation of the surface tension is much trickier because 
small crystallites are far from spherical and fluctuate considerably in shape for a given size. 
Crystallites become mostly spherical only well beyond the over-critical size. 

Figure U\ compares the free energy for these two methods: from the equilibrium proba- 
bility (Eq. |2J) and from CNT predictions (Eq. EJ). Following standard practice, the surface 
energy parameter a is fitted in order to obtain the best agreement with the first method. 
The Helmholtz free energy difference between the crystalline and liquid phases is computed 
as follows. 

The Gibbs free energy difference AG s i between solid and liquid states in NPT conditions 
at zero pressure is given by the difference in chemical potential Ap between the two phases. 
This quantity was computed by Broughton and Li^ 3 - and was found to be —7.697- 10~ 2 eV/at. 
However, we need the Helmholtz free energy difference AF s i at fixed density, which we can 
obtained by thermodynamical integration from the zero pressure results. Starting with the 
relation for the internal pressure [J^J — ~P, we use a thermodynamic integration for 
each phase (/- and s-Si) : 

AF= dF = - P(V) dV (4) 

JF! JV! 

The free energy difference between our system at zero pressure and at fixed density is 
computed by a Gaussian integration with five values: 

r PQO dV = (^^) E WiPQQ (5) 
Jvi \ A / , =1 

Vl =(^\«+(*±*\ (6) 



2 J V 2 

where Xi are the values for the Gaussian integration with their relative weight u>j. The initial 
volume, at zero pressure, for the liquid is Vxj = 18827.9 A 3 (2.467 g/cm 3 ) and the solid, 
V ltS = 20277.6 A 3 (2.29 g/cm 3 ); the final volume is V 2 = 20023.4 A 3 . Each point in the 
integral is simulated in NVT conditions at 1250 K for the liquid and solid. We equilibrate 
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our 1000-atoms system for 58 ps and then compute the mean pressure during 345 ps of 
simulation time. 

After integrating, we find a free energy difference per atom between the fixed density (p 
= 2.32 g/cm 3 ) and zero pressure system—, for liquid and solid state : 

AF S = 9.588 • 10" 4 eV/at. (7) 
AFi = 5.573 • 1(T 3 eV/at. (8) 

This gives a free energy difference between the liquid and solid phase at 2.32 g/cm 3 and 
1250 K of 

AF sl = AG sl (P = 0) + AF S - A F t (9) 
AF s i = -8.158 • 1(T 2 eV/at. (10) 

The const ant- volume correction is therefore only 6% of the zero pressure result of Broughton 
and Li. 

As can be seen in Fig. [7[ however, the CNT curve does not match the free energy data 
coming from P eq {n) in simulations. In order to find a better fit, the free energy difference 
AF s j between the solid and liquid state should be nine times lower than the value computed 
with the thermodynamic integration. 

We can verify the impact due to the choice of the order parameter on the free energy 
curve by re-analyzing the data using the criteria of Uttormark et al. ~. The resulting curve 
is also plotted in Fig. [7| and shows an even flatter curve, away from CNT results. We also 
repeated the simulation at 1250 K in the NPT ensemble at zero pressure and over 10 ns. In 
this situation, the trajectories do not crystallize — the largest crystallite reaches about 100 
atoms, well below the estimated critical size. The free energy distribution obtained from 
the cluster size distribution, while more curved than that for the NVT conditions, is still far 
from the CNT predictions (AG s i is about 5 times too low). 

The discrepancy between the two approaches clearly indicates that the classical nucleation 
theory does not fully capture the nucleation process in SW Si. We identify two possible 
sources of discrepancy. (1) As was demonstrated by Sastry and Angell recentl y 1 ^ 1 ? , SW Si 
undergoes a high-density to low-density liquid-liquid phase transition. The low density phase 
could be stabilized at higher temperature by the presence of a crystallite. In this case, it 
would be necessary to take into account two interfaces instead of one in the CNT equations. 
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FIG. 8: Free energy (divided by fc^T) of crystallites as a function of their size in NPT conditions 
(NPT) compared to NVT conditions (NVT). The simulation data are computed from the equilib- 
rium probability of presence for clusters with the basic blocks analysis. The CNT curve computed 
with the AG si of Broughton et Li (Fit A) is closer to the free energy data originating from P eq {n). 
However, a better fit (Fit B) requires a value five times lower. 

(2) The CNT fails because the critical nucleus is too small breaking the approximation 
of spherical crystallites. At this moment, we could not verify or disprove either of these 
possibilities. 

3. Lifetime of crystallites 

Beyond the free energy curve, we also analyze the dynamics of the crystallites present in 
the supercooled liquid. 

The lifetime probability of crystallites can be derived by following the kinetic approach 
of Zeldovich — . This approach predicts that the evolution of the clusters can be described 
by a diffusion equation of the form: 



where c(n,t) is the concentration of crystallites of size n at time t, D(n) the diffusion and 
&b T the boltzmann factor times the temperature. 

Van Kampen^ resolved the differential equation for small times by assuming the diffusion 
to be constant. Further approximating the potential as linear with respect to the cluster's 
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size, van Beijeren 3 * succeeded in finding a solution for longer times. This latter equation, 
which gives the distribution function for first arrival at size rif, for crystallites starting from 
size no is a well-known results 35 that confirms van Kampen short term behavior and contains 
an additional friction term e~ Uot which becomes important for longer times: 



)2 (AG(n,)-AG(» )) 

P (n, nf , t) = n ° e -^iD^e ^ e ~ Vot (12) 

V4:TCDt 3 

where 

This probability distribution is formally known as an inverse Gaussian distribution (or 
inverse normal, Wald). It was first derived independently by Schrodinger^ and Smolu- 
chowski 37 to describe Brownian motion in systems with a drift velocity. Hence, the devel- 
opment of a crystallite can be represented as a random walk in a field of force AF through 
different size classes where small clusters have a strong tendency to dissolve into the liquid 
(a drift to n/ — > 0) and super-critical nano-crystals tend to growth to macroscopic size ( 
rif — > oo). Since we do not have all the information on the free energy of crystallites AG(n) 
(see section UlI B 2j) and the diffusion constant, it is not possible to use directly Eq. IT2*1 to 
compare the lifetime behavior of clusters in the supercooled liquid during nucleation. How- 
ever, we can circumvent the difficulty by writing the inverse Gaussian distribution under a 
parametric form where A represents the mean and A/B 3 , the variance: 



= -(L-np) k B T 

D (AG(L) - AG(n )) 1 ; 



A 3 -2(L - n ) ( k B T 



(16) 



B D 2 \AG(L) - AG(n ). 

We compute the mean lifetime and variance for crystallites reaching the same maximum 
size in order to determine the theoretical distribution and compare with the lifetime prob- 
abilities from numerical simulations. Because large clusters are not encountered frequently, 
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FIG. 9: Lifetime distribution of clusters reaching a maximum size of 10 and 30 atoms. Comparison 
between the inverse gaussian distribution (theory) and the simulations data (simulation) with an 
uncertainty of ± 0.58 ps. 

the amount of data collected over all MD simulations remains small for the lilfetime of 
cluster near the critical size. In Fig. El the lifetime distributions determined by the Eq. El 
and the simulations data are in good agreement for small crystallites ensuring that cluster 
nucleation is well described by the inverse Gaussian distribution. The mean lifetimes for 
crystallites reaching an under-critical size of 10 or 30 atoms is 1.32± 0.6 and 3.79 ± 0.6 ps 
respectively, with a variance of 1.73 and 10.01, although some rare clusters last until 30 and 
50 ps (not shown). As would be expected, the mean lifetime increase with the size. 

Although some approximations have been made to obtain the lifetime probabilities of 
clusters by the inverse gaussian distribution and from the simulations data, the results 
are conclusive for crystallites reaching relatively small size. Since small clusters developed 
themselves in a confined range of size, we believe that the free energy difference can be 
approximated by a linear relation to the crystallite size and the diffusion kept constant. 

IV. CONCLUSIONS 

There has been a lot of interest recently regarding the nature of the liquid-solid transition 
in oriented liquid such as water and tetrahedral semiconductors. In many systems, it appears 
that there exists a high-density to low-density liquid transition often leading to a glassy or 
amorphous phas o 1 ^ 1 ^?? . Here, we reported results on a study of nucleation in liquid Si above 
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the HDL to LDL transition. 

We find that homogeneous nucleation takes place on a time scale of about 10 ns in a 
large enough system at constant volume. Using a topological order parameter, it is possible 
to follow the evolution of the crystallites through the crystallization process. Based on this 
analysis, we estimate the critical size to be around 175 atoms, within the limits of previous 
estimation of Uttormark et al. Surprisingly, the critical cluster, the one that will eventually 
crystallize the whole system, can survive at under-critical size for a long time (up to 900 ps 
or more) before it starts to grow steadily. Although the cluster's center of mass does not 
move significantly, there is a fluctuation in the composition of the cluster, as atoms move 
from the liquid to the crystallite and vice-versa, even once the crystallite has reached an 
over-critical size. 

A comparison of the simulation results with the classical nucleation theory indicates that 
the general behavior of the nucleation process is in agreement with CNT. For example, we 
find that the lifetime distribution of clusters reaching a specific maximum size follows the 
inverse gaussian distribution predicted recently^, supporting the description of the cluster 
growth as a random walk in the presence of a force field associated with the free energy. 
However, the details of the nucleation free energy differ significantly from the theoretical 
predictions. While the specific origin of this discrepancy remains open, we suggest that it 
could be caused by the presence of a low-density liquid at the interface between the crystal 
and the normal liquid or by the small size of the critical nucleus. More studies are required 
to fully address this problem. 
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